# 读取文件,表格行列转置,index变成cell对特征进行pca
file=open("romanov.rds",mode='r',encoding='cp936')
import pandas as pd
data=pd.read_csv(file,sep='\t')
data_T=pd.DataFrame(data.T,columns=data.index,index=data.columns)

# pca降维
from sklearn.decomposition import PCA
pca=PCA(n_components=16)
pca.fit(data_T)
data_Tpca=pca.transform(data_T)
# print(pca.explained_variance_)
# print(pca.explained_variance_ratio_)

# Kmeans聚类
# from sklearn.cluster import KMeans
# kmeans=KMeans(n_clusters=7,random_state=0)
# kmeans.fit(data_Tpca)
# lable=kmeans.labels_

# 层次聚类,不同簇之间的距离使用欧几里得距离,簇中点间的平均距离作为簇距离average/complete/single
# from sklearn.cluster import AgglomerativeClustering
# clu=AgglomerativeClustering(n_clusters=7,affinity='euclidean',linkage='complete')
# clustering=clu.fit(data_Tpca) # x 为array 样本特征
# lable=clustering.labels_  # 输出每个点的label

# DBSCAN
from sklearn.cluster import DBSCAN
dbscan=DBSCAN(eps=12.8,min_samples=3,metric='euclidean')
dbscan.fit(data_Tpca)
lable=dbscan.labels_
# dbscan.predict(lable)
print(lable)

# 保存聚类结果,用来导出结果表格
lable_t=lable.reshape(-1,1)
lable_e=pd.DataFrame(lable_t,index=data_T.index,columns=['lable'])
# lable_e.to_csv('klables.csv')

# 计算ARI和NMI
tdata=pd.read_csv('romanov.rds_label',sep='\t')
tdata=tdata.drop(columns='cell_type2')
tlable=tdata.to_numpy()
tlable=tlable.flatten()
# print(tlable)
from sklearn import metrics
print(metrics.adjusted_rand_score(tlable,lable))
NMI=metrics.normalized_mutual_info_score(tlable,lable)
print(NMI)

# 可视化聚类结果（选取pca降维结果中特征方差值最大的两列作为横纵坐标）
import matplotlib.pyplot as plt
red_x,red_y=[],[]
blue_x,blue_y=[],[]
green_x,green_y=[],[]
black_x,black_y=[],[]
cyan_x,cyan_y=[],[]
redx_x,redx_y=[],[]
yellow_x,yellow_y=[],[]


# print(reduced_x)
for i in range(len(data_Tpca)):
  if lable[i]==0:
    red_x.append(data_Tpca[i][0])
    red_y.append(data_Tpca[i][1])
  elif lable[i]==1:
    blue_x.append(data_Tpca[i][0])
    blue_y.append(data_Tpca[i][1])
  elif lable[i]==2:
    black_x.append(data_Tpca[i][0])
    black_y.append(data_Tpca[i][1])
  elif lable[i]==3:
    green_x.append(data_Tpca[i][0])
    green_y.append(data_Tpca[i][1])
  elif lable[i]==4:
    cyan_x.append(data_Tpca[i][0])
    cyan_y.append(data_Tpca[i][1])
  elif lable[i]==-1:
    redx_x.append(data_Tpca[i][0])
    redx_y.append(data_Tpca[i][1])
  else:
    yellow_x.append(data_Tpca[i][0])
    yellow_y.append(data_Tpca[i][1])
plt.scatter(red_x,red_y,c='m',marker='.')
plt.scatter(blue_x,blue_y,c='b',marker='.')
plt.scatter(green_x,green_y,c='g',marker='.')
plt.scatter(black_x,black_y,c='k',marker='.')
plt.scatter(cyan_x,cyan_y,c='c',marker='.')
plt.scatter(redx_x,redx_y,c='r',marker='x')
plt.scatter(yellow_x,yellow_y,c='y',marker='.')
plt.show()

# print(tdata)
# 可视化真实lable的结果
# import matplotlib.pyplot as plt
# red_x,red_y=[],[]
# blue_x,blue_y =[],[]
# green_x,green_y =[],[]
# black_x,black_y=[],[]
# cyan_x,cyan_y=[],[]
# magenta_x,magenta_y=[],[]
# yellow_x,yellow_y=[],[]
#
#
# for i in range(len(data_Tpca)):
#   if tlable[i]=='oligos':
#     red_x.append(data_Tpca[i][0])
#     red_y.append(data_Tpca[i][1])
#   elif tlable[i]=='astrocytes':
#     blue_x.append(data_Tpca[i][0])
#     blue_y.append(data_Tpca[i][1])
#   elif tlable[i]=='ependymal':
#     black_x.append(data_Tpca[i][0])
#     black_y.append(data_Tpca[i][1])
#   elif tlable[i]=='microglia':
#     green_x.append(data_Tpca[i][0])
#     green_y.append(data_Tpca[i][1])
#   elif tlable[i]=='vsm':
#     cyan_x.append(data_Tpca[i][0])
#     cyan_y.append(data_Tpca[i][1])
#   elif tlable[i]=='endothelial':
#     magenta_x.append(data_Tpca[i][0])
#     magenta_y.append(data_Tpca[i][1])
#   else:
#     yellow_x.append(data_Tpca[i][0])
#     yellow_y.append(data_Tpca[i][1])
# plt.scatter(red_x,red_y,c='r',marker='.')
# plt.scatter(blue_x,blue_y,c='b',marker='.')
# plt.scatter(green_x,green_y,c='g',marker='.')
# plt.scatter(black_x,black_y,c='k',marker='.')
# plt.scatter(cyan_x,cyan_y,c='c',marker='.')
# plt.scatter(magenta_x,magenta_y,c='m',marker='.')
# plt.scatter(yellow_x,yellow_y,c='y',marker='.')
# plt.show()